import numpy as np
import cv2
import matplotlib.pyplot as plt
from scipy.ndimage import convolve as convolveim
import scipy.ndimage.filters as filters
import scipy.ndimage as ndimage
def gauss(size,sigma):
"""
This function will generate 3 filters given the size of the filter and sigma of Gaussian:
1: gaussian filter;
2: derivative of gaussian filters in x and y direction.
"""
# define the x range
x_ax = np.arange(0,size) - size/2 + 0.5
################################################ TODO ###############################################
# make 1D gaussian filter
gauss = np.exp(-x_ax**2/(2*sigma**2))
# Compose 2D Gaussian filter from 1D, using the separability property of 2D Gaussian
gauss2 = gauss.reshape(size,1) * gauss.reshape(1,size)
# Normalize the filter so that all coefficients sum to 1
gauss1 = gauss2/np.sum(gauss2)
# Create derivatives of gaussian
gauss1_dx = np.matrix(np.zeros((np.shape(gauss1)), dtype="float32"))
gauss1_dy = np.matrix(np.zeros((np.shape(gauss1)), dtype="float32"))
for j in range(0, len(x_ax)):
# derivative filter in x
gauss1_dx[:, j] = (gauss1[:, j] * (-x_ax[j])/(sigma*sigma)).reshape(size,1)
################################################ TODO ###############################################
# similarly define the difference in y
gauss1_dy[j, :] = (gauss1[j, :] * (-x_ax[j])/(sigma*sigma)).reshape(1,size)
return gauss1,gauss1_dx, gauss1_dy
# Visualize the filters you created to make sure you are working with the correct filters
gauss1, gauss1_dx, gauss1_dy = gauss(21,3)
plt.figure()
plt.subplot(1,3,1)
plt.imshow(gauss1)
plt.subplot(1,3,2)
plt.imshow(gauss1_dx, cmap='gray')
plt.subplot(1,3,3)
plt.imshow(gauss1_dy, cmap='gray')
plt.show()
# visualize the filters in 3D
from mpl_toolkits.mplot3d import Axes3D
x = np.arange(0,21)
y = np.arange(0,21)
X, Y = np.meshgrid(x, y)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(X, Y, gauss1, cmap='viridis')
plt.show()
def harris(Ix, Iy , input_image,N):
"""
The input to this function are the gradient images in x and y directions, the original image and N
The function will output two arrays/lists x and y which are the N points with largest harris values,
and an image of the harris values.
"""
l, m = np.shape(input_image)
################################################ TODO ###############################################
#Forming 3 images
#Ix square
Ix2 = Ix*Ix
#Iy square
Iy2 = Iy*Iy
#Ix*Iy
Ixy = Ix*Iy
# Smooth image Ix2, Iy2, Ixy with Gaussian filter with sigma=2, size=7.
# Get the gauss filter for smoothing (reuse what you have)
gauss_smooth, _, _ = gauss(7,2)
Ix2_smooth = convolveim(Ix2, gauss_smooth, mode='nearest')
################################################ TODO ###############################################
# CONVOLVE as shown above
Iy2_smooth = convolveim(Iy2, gauss_smooth, mode='nearest')
Ixy_smooth = convolveim(Ixy, gauss_smooth, mode='nearest')
# By doing this, Ix2_smooth, Iy2_smooth, Ixy_smooth are the three values needed to calculate
# the A matrix for each pixel.
################################################ TODO ###############################################
# Write code segment to find N harris points in the image
# Refer to the page 17 of slides on features for the equation
det = Ix2_smooth*Iy2_smooth - Ixy_smooth**2
trace = Ix2_smooth + Iy2_smooth
H = det - 0.06*trace**2
################################################ TODO ###############################################
# Save a copy of the original harris values before detecting local max
H0 = H
# Detect local maximum over 7x7 windows
local_max_win = 7
a = int(np.floor(local_max_win/2))
H = np.pad(H,((a,a),(a,a)), 'constant')
# Initialize a mask to be all ones. The mask corresponds to the local maximum in H
H_max = np.ones(H.shape)
for i in range(a,l+a):
for j in range(a,m+a):
# Take a WxW patch centered at point (i,j), check if the center point is larger than all other points
# in this patch. If it is NOT local max, set H_max[i,j] = 0
patch = H[i-a:i+a+1,j-a:j+a+1]
if H[i,j] < np.max(patch):
H_max[i,j] = 0
# Multiply the mask with H, points that are not local max will become zero
H = H_max*H
H = H[a:-a,a:-a]
# Find largest N points' coordinates
# Hint: use np.argsort() and np.unravel_index() to sort H and get the index in sorted order
# Sort the harris values in descending order
idx = np.argsort(H.flatten())[::-1][:N]
# Get the coordinates of the N largest harris values
x, y = np.unravel_index(idx, H.shape)
# x,y should be arrays/lists of x and y coordinates of the harris points.
return x,y,H0
##### IMPORTANT: Convert your image to float once you load the image. ######
input_image = cv2.imread('9.png',0).astype('float')
img = cv2.normalize(input_image, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
################################################ TODO ###############################################
# Generating the gaussian filter
sigma = 1
size = int(4*sigma + 1)
# Function call to gauss
gauss_filt, gauss_filt_dx, gauss_filt_dy = gauss(size,sigma)
################################################ TODO ###############################################
# Convolving the filter with the image
# Convolve image with dx filter
Ix = convolveim(input_image,gauss_filt_dx,mode ='nearest')
# Convolve image with dy filter
Iy = convolveim(input_image,gauss_filt_dy,mode ='nearest')
x,y,H0 = harris(Ix, Iy ,input_image,100)
################################################ TODO ###############################################
# Plot: Ix, Iy, Original image with harris point labeled in red dots, H0 harris value image
# Hint: you may use "plt.plot(y,x, 'ro')" # Note: x is vertical and y is horizontal in our above definition
# But when plotting the point, the definition is reversed
fig, axs = plt.subplots(2, 2, figsize=(12, 12))
axs[0, 0].imshow(Ix, cmap='gray')
axs[0, 0].set_title('Ix')
axs[0, 1].imshow(Iy, cmap='gray')
axs[0, 1].set_title('Iy')
axs[1, 0].imshow(input_image, cmap='gray')
axs[1, 0].plot(y, x,'ro')
axs[1, 0].set_title('Original image with Harris points')
axs[1, 1].imshow(H0, cmap='gray')
axs[1, 1].set_title('Harris values')
plt.show()
Write a program that can generate SIFT descriptor for each detected feature point using your program in Prob. 1. You may follow the following steps:
def histo(theta4,mag4):
"""
theta4: an array of quantized orientations, with values 0,1,2...7.
mag4: an array of the same size with magnitudes
"""
temp = np.zeros((1,8),dtype='float32')
################################################ TODO ###############################################
# write code segment to add the magnitudes of all vectors in same orientations
for i in range(8):
temp[0,i] = np.sum(mag4[theta4==i])
# temp should be a 1x8 vector, where each value corresponds to an orientation and
# contains the sum of all gradient magnitude, oriented in that orientation
return temp
def descriptor(theta16,mag16):
"""
Given a 16x16 patch of theta and magnitude, generate a (1x128) descriptor
"""
filt,_,_ = gauss(16,8)
mag16_filt = mag16*filt
# array to store the descriptor. Note that in the end descriptor should have size (1, 128)
desp = np.array([])
################################################ TODO ###############################################
# Make function call to histo, with arguments theta16 and mag16_fil
# This is used for find the location of maximum theta
histo16 = histo(theta16,mag16_filt)
maxloc_theta16 = np.argmax(histo16)
for i in range(0,16,4):
for j in range(0,16,4):
################################################ TODO ###############################################
# Use histo function to create histogram of oriantations on 4x4 pathces in the neighbourhood of the harris points
# You should shift your histogram for each cell so that the dominant orientation of the 16x16 patch becomes the first quantized orientation
# You should update the variable desp to store all the orientation magnitude sum for each sub region of size 4x4
histo4 = histo(theta16[i:i+4,j:j+4],mag16_filt[i:i+4,j:j+4])
desp = np.append(desp,np.roll(histo4,-maxloc_theta16))
################################################ TODO ###############################################
# normalize descriptor, clip descriptor, normalize descriptor again
desp = desp/np.linalg.norm(desp)
desp = np.clip(desp,0,0.2)
desp = desp/np.linalg.norm(desp)
desp = np.matrix(desp)
return desp
def part_B(input_image):
# Normalize the image
img = cv2.normalize(input_image, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
# Generate derivative of Gaussian filters, using sigma=1, filter window size=4*sigma+1
sigma = 1
_, filt_dx, filt_dy = gauss(int(4*sigma + 1), sigma)
################################################ TODO ###############################################
# Image convolved with filt_dx and filt_dy
img_x = cv2.filter2D(img, -1, filt_dx)
img_y = cv2.filter2D(img, -1, filt_dy)
# Calculate magnitude and theta, then quantize theta.
mag = np.sqrt(img_x ** 2 + img_y ** 2)
theta = np.arctan2(img_y, img_x)
theta = (theta/(2*np.pi))*360
theta = theta*(theta>=0) + (360+theta)*(theta < 0)
################################################ TODO ###############################################
# Quantize theta to 0,1,2,... 7, see instructions above
q = 45
N = 8
theta_q = np.floor(theta/q)
theta_q[theta_q==N] = 0
theta_q = theta_q.astype('uint8')
################################################ TODO ###############################################
# Call harris function to find 100 feature points
x,y,_ = harris(img_x, img_y ,input_image,100)
# Pad 15 rows and columns. You will need this extra border to get a patch centered at the feature point
# when the feature points lie on the original border of the image.
theta_q = cv2.copyMakeBorder(theta_q.astype('uint8'), 7,8,7,8, cv2.BORDER_REFLECT)
################################################ TODO ###############################################
mag = cv2.copyMakeBorder(mag, 7,8,7,8, cv2.BORDER_REFLECT) # similarly add border to the magnitude image
final_descriptor = np.zeros((1,128))
final_points = np.vstack((x,y)).T
for i in range(final_points.shape[0]):
# Since you have already added 15 rows and columns, now the new coordinates of the feature points are (x+8, y+8).
# Then the patch should be [x[i]:x[i]+16,y[i]:y[i]16]
# Your patch should be centered at the feature point.
theta_temp = theta_q[x[i]:x[i]+16,y[i]:y[i]+16]
# similarly, take a 16x16 patch of mag around the point
mag_temp = mag[x[i]:x[i]+16,y[i]:y[i]+16]
# function call to descriptors
temp2 = descriptor(theta_temp,mag_temp)
final_descriptor = np.vstack((final_descriptor,temp2))
# Initially, final descriptor has a row of zeros. We are deleting that extra row here.
final_descriptor = np.delete(final_descriptor,0,0)
final_descriptor = np.nan_to_num(final_descriptor)
final_descriptor = np.array(final_descriptor)
# Combine x,y to form an array of size (Npoints,2) each row correspond to (x,y)
# You could use np.hstack() or np.vstack()
final_points = np.hstack((x.reshape(-1,1),y.reshape(-1,1)))
return final_descriptor,final_points
input_image = cv2.imread('9.png',0).astype('float') # input image
# Visualization the results. Plot the feature point similiar to Part1 and plot SIFT features as bar
final_descriptor , final_points = part_B(input_image)
for i in range(0,10):
f, (ax1, ax2) = plt.subplots(1, 2)
ax1.imshow(input_image,cmap='gray')
ax1.autoscale(False)
ax1.plot(final_points[i][1],final_points[i][0], 'ro')
ax2.bar(np.arange(1,129),final_descriptor[i,:])
plt.show()
Finding corresponding points in two images based on SIFT descriptors.
img1 = cv2.imread('left.jpg',0).astype('float') # read left image
img_1 = cv2.normalize(img1, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
descriptor_1, keypoints_1 = part_B(img_1)
img2 = cv2.imread('right.jpg',0).astype('float') # read right image
img_2 = cv2.normalize(img2, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_32F)
descriptor_2, keypoints_2 = part_B(img_2)
# Display detected points in the two images
plt.figure(figsize = (10,10))
plt.imshow(img_1,'gray')
plt.plot(keypoints_1[:,1],keypoints_1[:,0],'ro',ms=3)
plt.figure(figsize = (10,10))
plt.imshow(img_2,'gray')
plt.plot(keypoints_2[:,1],keypoints_2[:,0],'ro',ms=3)
[<matplotlib.lines.Line2D at 0x22b3b3347f0>]
# write function to find corresponding points in image
def points_matching(kp1, descriptor1, kp2, descriptor2, threshold):
matched_loca = list() # list of all corresponding points pairs. Point pairs can be stored as tuples
################################################ TODO ###############################################
# Find matching points between img1 and img2 using the algorithm described in the above
# For distance measuring, you may use np.linalg.norm()
# You could implement it as nested loop for simplicity.
for i in range(descriptor1.shape[0]):
distance = np.linalg.norm(descriptor1[i] - descriptor2, axis=1)
if np.min(distance) < threshold:
matched_loca.append((kp1[i], kp2[np.argmin(distance)]))
return matched_loca
# Test different thresholds for the matching
for r in [0.95, 0.8, 0.775, 0.75, 0.725, 0.7, 0.675, 0.65, 0.5]:
matched_loca = points_matching(keypoints_1, descriptor_1, keypoints_2, descriptor_2, r)
final_image = np.concatenate((img_1,img_2),axis=1)
print('threshold: ', r)
print('number of corresponding poitnts found:',len(matched_loca))
################################################ TODO ###############################################
# Write code segment to draw lines joining corresponding points
# Use cv2.line() to draw the line on final_image
# Remember the x,y coordinate in numpy and OpenCV is opposite and you need to add image width for pt2
for pt in matched_loca:
pt1 = (int(pt[0][1]), int(pt[0][0]))
pt2 = (int(pt[1][1] + img_1.shape[1]), int(pt[1][0]))
cv2.line(final_image, pt1, pt2, (255, 255, 255), 1)
plt.figure(figsize=(15,15))
plt.imshow(final_image,cmap='gray')
plt.show()
threshold: 0.95 number of corresponding poitnts found: 85
threshold: 0.8 number of corresponding poitnts found: 43
threshold: 0.775 number of corresponding poitnts found: 38
threshold: 0.75 number of corresponding poitnts found: 34
threshold: 0.725 number of corresponding poitnts found: 32
threshold: 0.7 number of corresponding poitnts found: 30
threshold: 0.675 number of corresponding poitnts found: 28
threshold: 0.65 number of corresponding poitnts found: 26
threshold: 0.5 number of corresponding poitnts found: 9
Stitch two images into a panorama using SIFT feature detector and descriptor. In this part, you may use functions from cv2.
- sift = cv2.SIFT_create()
- skp = sift.detect(img,None)
- (skp, features) = sift.compute(img, skp)
In your report, show the left and right images, the left and right images with SIFT points indicated, the image that illustrates the matching line between corresponding points, the transformed left image, and finally the stitched image.
def drawMatches(imageA, imageB, kpsA, kpsB, matches, status, lw=1):
# initialize the output visualization image
(hA, wA) = imageA.shape[:2]
(hB, wB) = imageB.shape[:2]
vis = np.zeros((max(hA, hB), wA + wB), dtype="uint8")
vis[0:hA, 0:wA] = imageA
vis[0:hB, wA:] = imageB
# loop over the matches
for ((trainIdx, queryIdx), s) in zip(matches, status):
# only process the match if the keypoint was successfully
# matched
if s == 1:
# draw the match
ptA = (int(kpsA[queryIdx][0]), int(kpsA[queryIdx][1]))
ptB = (int(kpsB[trainIdx][0]) + wA, int(kpsB[trainIdx][1]))
cv2.line(vis, ptA, ptB, (255, 255, 255), lw)
# return the visualization
return vis
img1 = cv2.imread('left.jpg',0) # read left image
img2 = cv2.imread('right.jpg',0) # read right iamge
# Depending on your OpenCV version, you could set up SIFT differently
sift = cv2.SIFT_create()
# sift = cv2.xfeatures2d.SIFT_create()
################################################ TODO ###############################################
# Use sift.detect to detect features in the images
kp1 = sift.detect(img1,None)
kp2 = sift.detect(img2,None)
# Visualize the keypoints
img1_kps = cv2.drawKeypoints(img1,kp1,None,flags = cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
img2_kps = cv2.drawKeypoints(img2,kp2,None,flags = cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS)
plt.figure(figsize=(15,15))
plt.subplot(121)
plt.imshow(img1_kps)
plt.title('image 1 with keypoints')
plt.subplot(122)
plt.imshow(img2_kps)
plt.title('image 2 with keypoints')
plt.show()
################################################ TODO ###############################################
# Use sift.compute to generate sift descriptors/features
(kp1, features1) = sift.compute(img1,kp1)
(kp2, features2) = sift.compute(img2,kp2)
kp1 = np.float32([kp.pt for kp in kp1])
kp2 = np.float32([kp.pt for kp in kp2])
matcher = cv2.DescriptorMatcher_create("BruteForce")
################################################ TODO ###############################################
# Use knnMatch function in matcher to find corresonding features
# For robustness of the matching results, we'd like to find 2 best matches (i.e. k=2 for knnMatch)
# and return their matching distances
rawMatches = matcher.knnMatch(features1, features2, 2)
matches = []
# Now we validate if the matching is reliable by checking if the best maching distance is less than
# the second matching by a threshold, for example, 20% of the 2nd best maching distance
for m in rawMatches:
################################################ TODO ###############################################
# Ensure the distance is within a certain ratio of each other (i.e. Lowe's ratio test)
# Test the distance between points. use m[0].distance and m[1].distance
if len(m) == 2 and m[0].distance < m[1].distance * 0.65:
matches.append((m[0].trainIdx, m[0].queryIdx))
ptsA = np.float32([kp1[i] for (_,i) in matches])
ptsB = np.float32([kp2[i] for (i,_) in matches])
################################################ TODO ###############################################
### Similar to what we did in part C
### Create an image img_match that shows the matching results by drawing lines between corresponding points.
img_match = drawMatches(img1, img2, kp1, kp2, matches, status=)
for p1, p2 in zip(ptsA, ptsB):
cv2.line(img_match, (int(p1[0]), int(p1[1])), (int(p2[0]), int(p2[1])), (255, 255, 255), 1)
plt.figure(figsize=(15,15))
plt.imshow(img_match,'gray')
plt.title('Matching points (Before RANSAC)')
plt.show()
################################################ TODO ###############################################
# Find homography with RANSAC
(H, status) = cv2.findHomography(...)
img_ransac = drawMatches(img1,img2,kp1,kp2,matches,status)
plt.figure(figsize=(15,15))
plt.imshow(img_ransac,'gray')
plt.title('Matching points (After RANSAC)')
plt.show()
################################################ TODO ###############################################
# fill in the arguments to warp the second image to fit the first image.
# For the size of the resulting image, you can set the height to be the same as the original, width to be twice the original.
# First transform the right image, then fill in the left part with the orignal left image
result = cv2.warpPerspective(...)
# For blending, you could just overlay img1 to the corresponding positions on warped img2
result[...] = img1
plt.figure(figsize=(20,40))
plt.imshow(result,'gray')
plt.title('Stitched image')
--------------------------------------------------------------------------- error Traceback (most recent call last) Cell In[14], line 3 1 ################################################ TODO ############################################### 2 # Find homography with RANSAC ----> 3 (H, status) = cv2.findHomography(...) 5 img_ransac = drawMatches(img1,img2,kp1,kp2,matches,status) 6 plt.figure(figsize=(15,15)) error: OpenCV(4.6.0) :-1: error: (-5:Bad argument) in function 'findHomography' > Overload resolution failed: > - findHomography() missing required argument 'dstPoints' (pos 2) > - findHomography() missing required argument 'dstPoints' (pos 2) > - findHomography() missing required argument 'dstPoints' (pos 2) > - findHomography() missing required argument 'dstPoints' (pos 2)
from __future__ import print_function
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
matched_idx = status.nonzero()[0]
def visualize_match(x):
"""
This function visualize the matches interactively
You could change the visualization of the matching keypoints by toggling a bar
Need to have the matches and status ready
matches: coarse matching results obtained from knnMatch
status: the refined matching results provided by cv2.findHomography,
the positive match determined by RANSAC is marked with 1,
while the negative match is marked with 0
"""
idx = matched_idx[x]
img_ransac = drawMatches(img1,img2,kp1,kp2,matches[idx:idx+1],status[idx:idx+1], lw=2)
plt.figure(figsize=(25,25))
plt.imshow(img_ransac,'gray')
plt.title('Matching points (After RANSAC)')
plt.show()
interact(visualize_match, x=widgets.IntSlider(min=0, max=len(matched_idx)-1, step=1, value=100));